/*
 * lib.h
 *
 *  Created on: 2010.06.26
 *      Author: Tadas
 */

#ifndef LIB_H_
#define LIB_H_

#include <cmath>

using namespace std;

#define CHANGE(a, b) {temp = (a); (a) = (b); (b) = temp;}
#define TINY 1.0e-20

double * allocate1d(int imax)
{
	double * mas;
	int i;

	mas = new double[imax];
	assert (mas != NULL);
	for (i=0 ; i<imax ; i++) mas[i] = 0;

	return mas;
}

void deallocate2d(double **mas, int eil)
{
	for (int i=0; i<eil; i++)
		delete [] mas[i];

	delete [] mas;
}

void deallocate1d (double * mas){	
	delete [] mas;	
}

double **alloc(int eil, int stl)
{
	double **mas;

	mas = new double *[eil];
	assert(mas!=NULL);

	for (int i=0; i<eil; i++)
		{
			mas[i] = new double [stl];
			assert(mas[i]!=NULL);
		}

	return mas;
}

int gausjord(double **a, double *b, double *x, int n)
{
	int *index;
	int i, icol, irow, j, k, sign;
	double big, dum, temp;

	index = new int[n];

	for (i=0; i<n; i++) index[i]=i;

	for (i=0; i<n; i++)
	{
		big = 0.0;

		for (j=i; j<n; j++)
		{
			for (k=i; k<n; k++)
			{
				if (fabs(a[j][k]) >= big)
				{
					big = fabs(a[j][k]);
					sign =1;
					if (a[j][k] < 0) sign = -1;
					irow = j;
					icol = k;
				}
			}
		}
		if (big == 0.0)
		{
			printf("guasjord: Singuliarine matrica");
			return 1;
		}
		if (icol!=i)
		{
			CHANGE(index[i], index[icol]);
			for (j = 0; j<n; j++) CHANGE(a[j][i], a[j][icol]);
		}
		if (irow!=i)
		{
			CHANGE(b[i], b[irow]);
			for (j = i; j<n; j++) CHANGE(a[i][j], a[irow][j]);
		}
		for ( j = i; j<n; j++) a[i][j] /= big*sign;
		b[i] /= big * sign;

		for (j = i+1; j<n; j++)
		{
			dum = a[j][i];
			for (k = i; k<n; k++) a[j][k] -= dum * a[i][k];
			b[j] -= dum * b[i];
		}
	}

	for (i = n-1; i>=0; i--)
	{
		temp = b[i];

		for (j=i+1; j<n; j++) temp -= a[i][j]*b[j];
		b[i] = temp / a[i][i];
	}

	for (i=0; i<n; i++) x[index[i]] = b[i];

	delete [] index;
	return 0;
}

int lowup(double **mas, int n){
	int i, j, k;
	double sum;
	for(j = 0; j < n; j++){
		for(i = 0; i <= j; i++){
			sum = mas[i][j];
			for(k = 0; k < i; k++){
				sum -= mas[i][k]*mas[k][j];
			}
			mas[i][j] = sum;
		}
		
		for( i = j + 1; i < n; i++){
			sum = mas[i][j];
			
			for(k = 0; k < j; k++){
				sum -= mas[i][k]*mas[k][j];
			}
			mas[i][j] = sum/mas[j][j];
		}
	}
	return 0;
}

void solve(double **mas, double *b, double *x, int n){
	double *y, sum;
	int i, j;
	
	y = allocate1d(n);
	y[0] = b[0];
	
	for(i = 1; i < n; i++){
		sum = 0.0;
		for( j = 0; j < i; j++){
			sum += mas[i][j] * y[j];
		}
		y[i] = (b[i] - sum);
	}
	
	x[n - 1] = y[n - 1]/mas[n - 1][n - 1];
	
	for (i = n - 2; i >= 0; i--){
		sum = 0;
		for(j = i + 1; j < n; j++){
			sum += mas[i][j]*x[j];
		}
		x[i] = (y[i] - sum)/mas[i][i];
	}
	
	deallocate1d(y);
}

void tridag(double *a, double *b, double *c, double *x, double *r, unsigned long n){
	unsigned long j;
	double bet, *san;
	
	san = allocate1d(n);
	
	if(b[0] == 0.0){
		cout << "Klaida funkcijoje tridag" << endl;
		return;
	}
	
	x[0] = r[0]/(bet = b[0]);
	
	for(j = 1; j < n; j++){
		san[j] = c[j - 1]/bet;
		bet = b[j] - a[j]*san[j];
		
		if( bet == 0.0){
			cout << "Klaida funkcijoje tridag 2" << endl;
			return;
		}
		
		x[j] = (r[j] - a[j]*x[j-1])/bet;
	}
	
	for(j = n - 2; j <= 0; j--){
		x[j] -= san[j+1] * x[j + 1];
	}
	
	deallocate1d(san);
}

int lumcmp( double **a, int * indx, int n){
	int i, imax, j, k;
	double big, dum, sum, temp;
	double * vv;
	
	vv = allocate1d(n);
	
	for(i = 0; i < n; i++){
		big = 0.0;
		for( j = 0; j < n; j++){
			if((temp = fabs(a[i][j])) > big){
				big = temp;
			}
		}
		if( big == 0.0){
			cout << "Singuliari matrica funkcijoje lumcpm" << endl;
			return 1;
		}
		vv[i] = 1.0/big;
	}
	
	for(j = 0; j < n; i++){
		for(i = 0; i < j; i++){
			sum = a[i][j];
			for(k = 0; k < i; k++){
				sum -= a[i][k] * a[k][j];
			}
			a[i][j] = sum;
		}
		
		big = 0.0;
		
		for(i = j; i < n; i++){
			sum = a[i][j];
			for(k = 0; k < j; k++){
				sum -= a[i][k] * a[k][j];
			}
			a[i][j] = sum;
			
			if((dum = vv[i] * fabs(a[i][j])) >= big){
				big = dum;
				imax = i;
			}
		}
		
		if( j != imax){
			for( k = 0; k < n; k++){
				CHANGE(a[imax][k], a[j][k]);
			}
			CHANGE(vv[imax], vv[j]);
		}
		
		indx[j] = imax;
		
		if(a[j][j == 0.0]){
			a[j][j] = TINY;
		}
		if(j != n - 1){
			dum = 1.0/a[j][j];
			for(i = j + 1; i < n; i++){
				a[i][j] *= dum;
			}
			
		}
		
	}
	
	deallocate1d(vv);
	return 0;
}

#endif /* LIB_H_ */
